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Abstract 

A temporal finite element method based on a mixed form of the Hamiltonian weak 
principle is developed for dynamics and optimal control problems. The mixed form of 
Hamilton’s weak principle contains both displacements and momenta as primary variables 
that are expanded in terms of nodal values and simple polynomial shape functions. Un- 
like other forms of Hamilton’s principle, however, time derivatives of the momenta and 
displacements do not appear therein; instead, only the virtual momenta and virtual dis- 
placements are differentiated with respect to time. Based on the duality that is observed 
to exist between the mixed form of Hamilton’s weak principle and variational principles 
governing classical optimal control problems, a temporal finite element formulation of the 
latter can be developed in a rather straightforward manner. Several well-known problems 
in dynamics and optimal control are illustrated. The example dynamics problem involves a 
time-marching problem. As optimal control examples, elementary trajectory optimization 
problems are treated. 


Introduction 


This paper examines a finite element approach to addressing optimal control problems. 
Hamilton’s principle has traditionally been used in analytical mechanics as a method of 
obtaining the equations of motion for dynamical systems. Bailey [1] followed by several 
others (see [2 - 4] for example) obtained direct solutions to dynamics problems using a 
form of Hamilton’s principle known as the law of varying action thus opening the door for 
its use in computational mechanics. ' 

More recently, it has been shown that expression of Hamilton’s law as a weak form 
(commdtily referred to as Hamilton’s weak principle or HWP) provides a powerful alter- 
native to numerical solution of ordinary differential equations in the time domain [5. 6j. 
The accuracy of the time-marching procedure derived in [5, 6] is competitive with stan- 
dard ordinary differential equation solvers. However, in order to derive an unconditionally 
stable algorithm in [5] ,\.reduced / selective integration had to be used. Further computa- 
tional advantages may be obtained in so-called mixed formulations of HWP in which the 
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generalized coordinates and momenta appear as independent unknowns [7]. Therein, an 
unconditionally stable algorithm emerges for the linear oscillator with exact integration. 
IiWP also has shown to be an ideal tool for obtaining periodic solutions for autonomous 
systems, as well as finding the corresponding transition matrix for perturbations about the 
periodic solution [8], These are complex two-point boundary value problems; its utility for 
these problems and its superior performance in mixed form strongly suggest that it could 
be used in optimal control problems. 

In this paper we show that HWP in mixed form provides a useful parallel to optimal 
control problems. Although finite elements have been applied to optimal control problems 
[9], the present formulation is believed to offer advantages over existing formulations. For 
example, it allows simple choices for shape functions in the finite element formulation [7] 
which simplifies element quadrature. 

We begin by summarizing HWP in both displacement and mixed forms. The mixed 
form is then applied to initial value problems in dynamics. Then, making use of the well- 
known analogy between dynamics and optimal control, we develop our Hamiltonian weak 
form for optimal control problems. Finally, we apply it to two relatively simple optimal 
control problems to illustrate its power. 

Hamilton’s Weak Principle for Dynamics 

In [1 - 8], one can see the potential of obtaining a direct solution in the time domain 
very much analogous to obtaining the solution of a beam deflection problem with the 
beam axial coordinate broken into several segments or finite elements. In the present 
case, however, it is the time interval which is broken into segments; thus, the term “finite 
elements in time” has been adopted by several investigators. 

Only recently has a mixed formulation of HWP been investigated as a computational 
tool for finite elements in time [7]. In this section we will formulate the mixed form of 
HWP and illustrate its application to dynamics problems. 


General Development 

To this aim, let us consider an arbitrary holonomic mechanical system. The configu- 
ration is completely defined by a set of generalized coordinates q. Further, let us denote 
with L(q,q,t ) the Lagrangean of the system, Q the set of nonconservative generalized 
forces applied to the system, and p = dL/dq the set of generalized momenta. Then the 
following variational equation, known as HWP [5], describes the real motion of the system 
between the two known times to and tf. 



8L dt -f- f 8q^ Q dt 

Jt 0 



i s 


to 


( 1 ) 


This particular variational equation is said to be in displacement form because it only in- 
volves the variation of q. Although this formulation has been shown to be of practical use in 
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dynamics [5, 6], an even more useful formulation may be derived if independent variations 
in both displacements and momenta are allowed, resulting in a mixed formulation. 

To derive the mixed formulation, we begin by defining the Hamiltonian as 


H(q,p,t) =p T q - L(q,q,t) 


( 2 ) 


Taking the variation of Eq. (2) and substituting for 8L in Eq. (1) results in 



q + 8q T p — 8H + 8q T Q) dt = Sq T p 


l i 


*0 


( 3 ) 


Integrating the first term of the integrand by parts yields 



T p — Sp T q — 8H + Sq 1 Q ) dt = ( 8q T p — 8p T q) 


h 

*0 


( 4 ) 


This is called a mixed formulation because it contains independent variations of q and p. 
It is also a weak form in the sense that all boundary conditions are of the natural or weak 
type. 

There are two main advantages of the mixed formulation over the displacement for- 
mulation. The first advantage is that the mixed formulation generally provides a more 
accurate solution for a given level of computational effort than does the displacement for- 
mulation. The second advantage is that a simpler choice of shape functions is allowed. 
Note in Eq. (4) that time derivatives of 8q and 8p are present. However, no time deriva- 
tives of q and p exist. Therefore, it is possible to implement linear shape functions for Sq 
and Sp and constant shape functions for q and p within each element. 

Let us break the time interval from to to tf into N equally spaced elements. The 
nodal values of these elements are t{ for * = 1, . . . , N + 1 where f 0 = and tf — t^ +1 . We 
define a nondimensional elemental time r as 

( 5 ) 

The linear shape functions for the virtual coordinates and momenta are 


t-ti _t- ti 

t i+1 - ti At 


8q = 8qi(l - t) + 8q i+1 r 
8p = 8pi( 1 - r) + 8p i+1 r 


( 3 ) 


For the generalized coordinates and momenta 
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and 
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qi if r = 0; 
flt if 0 < r < 1; 
(]i + 1 if r = 1 


f ^ if r = 0; 

p = < 7>i if 0 < r < 1; 

l pi + x if t = 1 


(S) 


When these shape functions are substituted into Eq. (4), one can either generate an im- 
plicit.. time-marching procedure for nonlinear problems or apply standard finite element 
assembly procedures to solve periodic or two-point boundary value problems [S] . When 
this formulation is applied to the linear oscillator, a time-marching algorithm emerges that 
is unconditionally stable. Higher-order (so-called p- version) elements could be developed, 
and they would certainly be attractive for linear problems or for nonlinear problems with 
nonlinearities of low order. For nonlinear problems in general, use of the crude shape 
functions allowable with the mixed method would seem to be more efficient than use of 
higher-order shape functions in a p- version. The reason for this is that, with the exception 
of the term involving Q, all element quadrature can be done by inspection regardless of 
the order of the nonlinearities. Detailed comparison of these methods is beyond the scope 
of the present paper but is being undertaken by the first author at this time. 

Example 1: Nonlinear Initial-Value Problems 

Applying the shape functions of Eqs. (6) - (S) to Eq. (4) for an initial value problem, 
we obtain a recursive set of nonlinear algebraic equations of the form 


fj {qiJli+uPuPi+i) = 0 j = 1, 2, . . . ,n 


( 9 ) 


where n is four times the number of degrees of freedom of the system. Eq. (9) can be 
solved by a Newton-Raphson method yielding an implicit time-marching procedure. The 
key advantage of using finite elements and a weak variational approach over numerical 
integration is that the solution (for linear problems) is stable for all time steps. In other 
words, no matter how large a time step is used, a finite approximation of the solution will 
be obtained. This unconditional stability is obtained without ad hoc procedures such as 
selective or reduced integration which are necessary in displacement formulations. 

We also point out that 2 qi = <n + q l+ \ and 2 pi = p t + p i+1 . Thus, it is possible to cut 
the number of equations and unknowns in half. This can be very useful for a multi-degree 
of freedom problem. 
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Consider a simple pendulum composed of a lumped mass m and a weightless bar of 
length l (see Fig. 1). The single generalized coordinate q is the angular displacement of the 
bar from the vertical. Denoting the kinetic energy of the system with I\ and the potential 
energy with V , then we may define the following: 


K = -mi 2 q 2 
2 

V = ro<7^(l — cos q) 
L — K — V 
p - dL/dq = ml 2 q 
H — pq - L 


( 10 ) 


There are no nonconservative forces 0 applied to this system. 

Substituting t = ti + rAt , along with Eq. (10), and substituting the shape functions 
defined in Eqs. (6) - (8) into Eq. (4) we obtain 


At J | ~ <? i ) Pi ~ m 0^ sir W* [%(! “ r ) + %+i r l 

- 8q i+1 pi+i + f>Pi+i<li+i + - 6p t q t = 0 

Carrying out the integration, we obtain the four independent equations 



mg£At sin qi 

= 0 

Pi - Pi - 

2 

Pi -Pi + 1 - 

'mgPAt sin qi 

= 0 

2 


<1 

< 

= 0 

<li 

Ql 2 mt- 


'to 

<1 

1 

= 0 

Qi+1 

<h 2ml 2 


(ID 


( 12 ) 


There are six unknowns; however, for an initial-value problem, we will specify q, and p, 
and solve for the remaining unknowns. Thus, Eq. (12) is of the form of Eq. (9). 

For this simple pendulum example, we will nondimensionalize the variables as follows. 
If we define io 2 = g fi, then a dimensionless time step At may be defined as At = ^-A t. 
Also, instead of solving directly for p , we will solve for the dimensionless p/m(ru. 

We will start our pendulum at q\ = 60° and pi = 0.0. The equations will be solved 
for At = 0.4, 0.8, and 1.6. Graphs of the solutions are shown in Fig. 2 and Fig. 3 and 
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At = 0.4 gives acceptable results for both displacement and angular velocity. Also, note 
that even the large 1.6 time step yields a finite approximation of the exact solution. 

Weak Principle for Optimal Control 

A definite analogy exists between the mixed formulation of HWP in dynamics and 
the first variation of the performance index in optimal control theory. Specifically, there 
is a similarity between the generalized coordinates and generalized momenta in dynamics 
and the states and costates in optimal control theory. 

Since the mixed formulation has proven to be so valuable in dynamics, we will for- 
mulate a weak variational formulation of the performance index. When deriving this 
formulation, we will keep two things in mind. First, the resulting formulation must satisfy 
the Euler- Lagrange equations and boundary conditions that have already been established 
in optimal control theory [10]. Second, whatever terms are necessary will be added to 
the formulation to transfer all strong boundary conditions to natural or weak boundary 
conditions. 


General Development 

We start with a performance index taken from Eq. (2.S.4) of [10]. The first variation 
of the performance index will be taken in a standard manner, except that states, costates, 
and controls will have arbitrary variations. Rather than setting the first variation equal 
to zero, however, it will be set equal to an expression which contains the terms that are 
necessary to transform all boundary conditions to the natural or “weak” type. The final 
weak form is then obtained by integration of this equation by parts in such a way that no 
time derivatives of states or costates appear. 

Consider a system defined by a set of n states x and a set of m controls u. Furthermore, 
let the system be governed by a set of state equations of the form x = f(x, u, t). We may 
denote elements of the performance index, J, with an integrand L(r,ix,t) and a discrete 
function of the states and time at the final time In addition, any terminal 

constraints placed on the states may be placed in the set of q functions and 

adjoined to the performance index by a set of q discrete Lagrange multipliers v. Finally, we 
may adjoin the state equations to the performance index with a set of Lagrange multiplier 
functions A (t) which will be referred to as costates. This yields a performance index of the 
form 


J 


y 

Jlo 


( \ T i - L- \ T f)dt- <f> 


h 


-u T ip\ 


(13) 


Taking the first variation of J and setting it equal to an expression chosen so that al! 
boundary conditions are of the weak type, one obtains 
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6 J = J ' [6A T (i - /) + Si T A -<5L- <5/ t A] dt 
- 6x^A/ - 6i/ 7 V|^ 

= %(A r i)| t/ -&ejfAo-6A r (i-*) *' 


where 



( 14 ) 


(15) 


The right hand side of Eq. (14) contains terms necessary to form all of the proper boundary- 
conditions as natural ones. The quantities x and A are discrete values of the states and 
costates at the initial (with subscript 0) and final times (with subscript /). Depending on 
the problem, these values will either be specified or left as unknowns; they need not coincide 
with the values of the states or costates taken on within the elements at the beginning and 
at the end of the time interval. This is yet another indication of the “weakness” of the 
formulation. 

From Eq. (14), we may directly write down a weak formulation. However, as stated 
earlier, one of the requirements of the formulation is that it satisfies the Euler- Lagrange 
equations and boundary conditions. To show this, we will integrate the Sx T term in 
Eq. (14) by parts, expand the variation of L, substitute Eq. (15), and group terms yielding 



+ — A — 8xq (a 0 — A 0 ^ 

+ SXj (xf — x f) — 6 A<f (x 0 — x 0 ) = 0 


where a;o,A o,£/, and A/ represent the values of those functions at the initial and final 
times, respectively. The coefficients of 8\ T , 5x T , and 8u T in the integrand are the three 
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correct Euler-Lagrange equations, Eqs. 2.8.15 - 2.8.17 from [10]. There axe also six trailing 
terms in Eq. (16) from which the boundary conditions can be determined. The first four 
of these equations and the sixth equation correspond to the correct boundary conditions 
in [10]. Namely, the requirement for the coefficient of 8v T to vanish yields Eq. (2.8.21). 
The requirement for the coefficient of 8tj to vanish is equivalent to Eq. (2.8.20). The 
requirement for the coefficient of SxJ to vanish shows that the final value of A approaches 

A/ as given in Eq. (15), which corresponds to Eq. (2.8.19). If A 0 is chosen as zero, the 
requirement for the coefficient of SxJ to vanish requires the initial value of A to approach 
zero as the accuracy of the approximation scheme is increased (such as from adding more 
finite elements); on the other hand, the requirement for the coefficient of SXj to vanish 
requires the initial value of x to approach Xo, in accordance with Eq. (2.8.18). Finally, the 
requirement for the coefficient of 6Xj to vanish demands that the final value of x approach 
the discrete value x /; this has no counterpart in [10] since the elements of x/ are usually 
unknown. 

Having satisfied our requirement that none of the fundamental equations are altered, 
we may now derive our weak formulation from Eq. (14). Recall from the mixed formulation 
for dynamics that we do not want time derivatives of q or p to appear in the formulation. 
Similarly, we do not wish for the time derivatives of x and A to appear in the present 
weak formulation for optimal control. Therefore, we will integrate the x term by parts in 
Eq. (14). The l'esulting equation is 




8i T X - Sx T 


'to 


( d ^V + ( d -L 




dx 


X 


— <5A r x — 8X 1 f — Su 


T , 


d YX\m T 


du J ) 


dt 


-6t,[L + \Tf + ?£ + ^)\ -SvT* 

C r A r jH A * T A m ^ 

— ox j Xj Aq o X j* x f — SXq xq =0 




(17) 


This is the governing equation for the weak Hamiltonian method for optimal control prob- 
lems of the form specified. It will serve as the basis for the finite element discretization 
described below. It should be noted that normal^ one will encounter various types of in- 
equality constraints on the states and controls in problems that deal with optimal control. 
This aspect has not been treated yet and will be the subject of future research. 

As in dynamics, we may choose linear shape functions for 8x and S\. We may choose 
piecewise constant shape functions for x and A. Thus, we will be working with shape 
functions similar to those of Eqs. (6) - (8). In addition, note that the time derivatives of 
u and 6u do not appear in the formulation. Thus, we let 
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U — U{ 

Su — Stii 


( 18 ) 


Plugging in the shape functions described for x,A, and iz, substituting r for and 
carrying out the simple integration from 0 to 1, we obtain 


N 




i=l 


' 2 \dx ) i 1 2 vax ). 


+ SXj ( Xi - ^ fi 


+ 6x1 


1 + 1 


At 


fn r 


df y , _At/a£V 

1 2 \dx) i 1 2\dx) i 




At — 
2 


<{ At 


Vasy, 1 


-«/(£+ 








(19) 


+ dxf A 0 - 8\f x 0 - Sx J N+1 \f + <5A^ +1 x/ = 0 


This is the general algebraic form of our Hamiltonian weak form for optimal control 
problems of the form specified. Eq. (19) is a system of algebraic equations. In fact, for 
A r elements, there are 2n(N + 1) + mN + q + 1 equations and 2n(N + 2) + mN + q + 1 
unknowns. Therefore, 2n of the 4 n endpoint values for the states and costates (io- Aq. x /. 
and A/) must be specified. In general, Xq (the initial conditions) is known in accordance 
with physical constraints. Also, A / can be specified in terms of other unknowns with the 
use of Eq. (15). Now we have the same number of equations as unknowns. 

Normally Eq. (19) can be solved by writing an explicit Jacobian and using a Newton- 
Raphson solution procedure. For the example problems which follow, the iteration pro- 
cedure will convei'ge quickly for a small number of elements with a trivial initial guess. 
Then, the answers obtained for a small number of elements can be used to generate initial 
guesses for a higher number of elements. Thus, a large number of elements can be solved 
with a very efficient run-time on the computer. 

Unfortunately, for some highly nonlinear problems, trivial initial guesses may not be 
adequate. This problem may be overcome by noting that the Sx{ and 8x {. p j equations in 
Eq. (19) are n(N + 1) equations that happen to be linear in the n(N + 1) unknown costates. 
Now, we need n(N + 1) fewer initial guesses and we can solve for the costates in terms 
of the other unknowns. This is particularly useful since generating initial guesses for the 
costates can be difficult. This has proven to be a very useful way of obtaining answers for 
highly nonlinear equations. 
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Example 2: A Fixed-Final-Time Problem 

As the first optimal control problem, we will examine the transfer of a particle to a 
rectilinear path (see Fig. 4). This is an example taken from [10], article 2.4. We will let 
xi and x 2 denote the position of the particle at a given time and x 3 and x 4 denote the 
particle’s velocity at a given time. The thrust angle u is the control and the particle has 
mass m and a constant acceleration a. 

The state equations are defined as 

-0010 
0 0 0 1 
0 0 0 0 
.0 0 0 0 

The final time T is fixed and we would like to maximize the final horizontal component 
of velocity. Thus, L = 0 and (f> — x 3 (T). There are also two terminal constraints on the 
states. These are that the particle arrive with a fixed final height ( h ) and that the final 
vertical component of velocity be zero. We do not care what the final horizontal component 
of position is. These constraints can be stated analytically as 


x + 



(20) 


?/> = 



( 21 ) 


The initial conditions are x(0) = x 0 = |_ 0 0 0 0 J r . Finally, we will eliminate the 

unknown A/’s by writing it in terms of other unknowns. In accordance with Eq. (15) 


Af = 



( 22 ) 


With L — 0, and t / fixed, along with all zero initial conditions on the states, then the 
general formulation of Eq. (19) takes the following form. 


N 


EH 


*=1 


+ SxJ +l 


->-¥©> 


At - 

+ sx i ( - — /« 


c \ T I — At — 
0A i+1 xi + — -fi 


- 6u 


r 


-§)> 


-Su T ^\ t) 


+ <5x^Aq — 6xJj +l Xf + <5Ajy +1 x/ = 0 


(23) 
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This produces a system of nonlinear algebraic equations whose size depends on the 
number of elements N. 

These equations are solved by writing down an explicit Jacobian and using a Newton- 
Raphson algorithm. Trivial initial guesses (that are never changed regardless of input 
parameters) are used for N = 2. These results are then used to obtain the initial guesses 
for arbitrary N by simple interpolation. In all results obtained to date for this problem, 
no additional steps are necessary to obtain results as accurate as desired. 

Representative numerical results for all four states versus dimensionless time t/T are 
presented in Figs. 5-8. For this example, we have taken h = 100.0, T = 20.0, and 
4 h/aT 2 — 0.8S97. The results for 2, 4, and 8 elements are plotted against the exact 
solution available in [10]. It can easily be seen that A T = 8 gives acceptable results for all 
the states. 

In Fig. 9, the control angle u versus dimensionless time t/T is presented. Once again, 
the results are seen to be excellent for N — 8. 

Three of the four costates are constants for all time and this method yields two of 
these exactly. The third costate is very close to the exact answer. The fourth costate 
corresponding to the vertical component of velocity A 4 is shown in Fig. 10. The results 
compare nicely with the exact results. 

A plot of the relative error of the performance index J = ;r 3 (T) and the endpoint 
multiplier rq versus the number of elements is shown in Fig. 11. It is seen to be nearly a 
straight line on a log-log scale. The slope of the line is about —2 which indicates that the 
error varies inversely with the square of TV. 

Notice in Fig. 11 that there is a bend in the endpoint multiplier curve. It is a common 
characteristic of mixed formulations for an error curve not to be monotonically decreasing. 


Example 3: A Free-Final-Time Problem 

The second optimal control problem is similar to Example 2 except that now the 
final time is free and we would like to obtain a given horizontal component of velocity 
(U) in the minimum time (see problem 9, article 2.7 of [10]). Our algorithm from the 
preceding example is readily modified to fit this problem by noting the following changes. 
The performance index is now the final time T; so <j> = 0 and L == 1. Also, there is an 
additional endpoint constraint on the states; namely that x 3 = U. With these changes, we 
have 


A / = 



(24) 


and 
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x 2 — h 
xi-U 

£4 


( 25 ) 


</> = 


Along with these changes to our equations, we also pick up the additional St f equation. 
The new system of equations is solved in the same manner described previous!}'. Again, 
trivial initial guesses are satisfactory for N = 2, and these answer are used to obtain initial 
guesses for arbitrary N . 

Representative numerical results for all four states versus dimensionless time t/T are 
presented in Figs. 12-15 for a case with ah/TJ 2 = 0.75. The results for 2, 4, and 8 elements 
ai'e plotted against the exact solution available in [10]. It can easily be seen that N = 8 
gives acceptable results for all the states. 

The control angle u versus dimensionless time t/T is presented in Fig. 16. Once again, 
the results are seen to be excellent for N = 8. 

As with the fixed time problem, three of the four costates are constants. The costate 
results are all as accurate or better then the costate depicted in Fig. 10. 

A plot of the relative error of the performance index J — T versus the number of 
elements N is shown in Fig. 17. As before, it is seen to be nearly a straight line on a log- 
log scale. The slope of the line is about —2 which indicates that the error varies inversely 
with the square of N . 

It should be noted that the computer time on a Cyber 990 is only about 2 seconds 
for N — 2, N = 4, and N = 8 and 3 seconds for N — 16. Thus, the run tim.- is relatively 
insensitive to N . 


Conclusions and Future Work 

In this report, a mixed form of Hamilton’s Weak Principle has been stated for dy- 
namics problems. Finite elements in time were applied to this formulation and a simple 
initial- value problem has been used to demonstrate the principles involved. A temporal 
finite element method based on a mixed form of the Hamiltonian weak principle was then 
developed for optimal control problems from the dynamics principles. It has been demon- 
strated that the mixed form allows for a simple choice of shape functions and is essentially 
self-starting. Two simple optimal control problems have been examined and the results are 
seen to compare excellently with the exact solutions for even a very few elements. Overall, 
the method provides very accurate results for the problems investigated to date with only 
a few elements and for minimal computational effort. 

Future research will be done in applying this method to practical problems such as 
development of on-board trajectory optimization algorithms for launch vehicles [11]. Such 
applications will require the reliability, efficiency, and self-starting characteristics illus- 
trated in the present approach. However, more work needs to be done to solve the equa- 
tions as efficiently as possible. One area with great potential for making the method more 
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efficient is to take advantage of the sparsity of the Jacobian. The method will also have to 
be extended to allow for inequality constraints on the states and controls. 
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Fig. 1: Nomenclature for Example 1. This is a simple pendulum composed of a lump ma 
777 and a weightless bar of length C. 






Fig. 2. Angular displacement <j versus dimensionless time for three values of the time step 
At. These solutions are being compared to the exact elliptic integral solution. Note that 
even a large time step yields a stable solution. 



Dimensionless Momentum 



Dimensionless Time 


Fig. 3: Dimensionless momentum p/m£ 2 u> versus dimensionless time for three values of the 
time step At. These solutions are being compared to the exact elliptic integral solution. 
Note that even a large time step yields a stable solution. 




Fig. 4: Nomenclature for Example 2. Transfer of a particle to a rectilinear path with a 
fixed final time. The problem is to maximize the final horizontal component of velocity 
subject to specified terminal conditions on the states. 
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Fig. 5: Dimensionless horizontal position Xi/h versus t/T. Note that the final horizon! 
component of position is not specified. 
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Fig. 7: Dimensionless horizontal velocity x$T/h versus t/T. Note that the performan' 
index J — x 3 (T). 
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Fig. 8: Dimensionless vertical velocity x 4 T/h versus t/T. The final vertical component of 
velocity is constrained to be zero at the final time. 
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Fig. 9: Control angle u versus t/T . The control is chosen so as to maximize the performance 
index J = x 3 (T). 
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Fig. 10: Vertical velocity costate A 4 versus tfT . The results for this costate are the least 
accurate of all the costates. 
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Fig. 12: Dimensionless horizontal position x^/h versus t/T. Note that the final horizontal 
component of position is not specified. 
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Fig. 13: Dimensionless vertical position x 2 /h versus t/T. The final height is constrained 
to be h at the final time. 
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Fig. 14: Dimensionless horizontal velocity x 3 /U versus t/T. The final horizontal com- 
ponent of velocity is constrained to be U . The problem is to reach this velocity in the 
- irdnumum time T. 
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Fig. 15: Dimensionless vertical velocity x 4 /l 
velocity is constrained to be zero at the fina 
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J versus t/T. The final vertical component of 
1 time. 
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Fig. 16 : Control angle u versus t/T. The accuracy of the initial control u{t 0 ) determine 
the accuracy of the costates. 
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Fig. 17: Relative e 
with the square of . 


